Introducción a la Psicometría Bayesiana
February 27, 2025
{blavaan}
{cmdstanr}¿Cómo incluimos nuestra incertidumbre en un modelo estadístico?
Si creo que vuestra edad está entre 22 y 26 años, puedo expresar esta creencia usando una distribución de probabilidad. Por ejemplo, puedo utilizar la distribución normal para reflejar lo que creo:
\[ \Pr(\mu_{\text{edad}})\sim\mathcal{N}\left(\color{royalblue}{\mu = 24},\,\color{firebrick}{\sigma = 1}\right) \]
Esto significa que espero que la media de edad esté en torno a 24 años y que, al mismo tiempo, creo que la edad del 95% de vosotros estará, aproxiadamente, entre 22 y 26 años:
\[ \Pr\left(\mu_{\text{edad}} \mid D\right) \propto \Pr\left(D \mid \mu_{\text{edad}}\right) \times \Pr(\mu_{\text{edad}}) \]
Pierre-Simon Laplace
La inferencia Bayesiana es sentido común expresado en números.
Dennis Lindley
La única estadística buena.
Frecuentismo
Bayesiano
¿Cuál de las dos es más antigua?
Formalmente, la estadística bayesiana, que surge en el siglo XVIII.
Thomas Bayes
El único procedimiento estadístico que es coherente, es decir, que evita afirmaciones que sean internamente inconsistentes.
El hombre de la foto no es Thomas Bayes
Aun así, este autorretrato sigue usándose en los cursos de estadística Bayesiana.
Cada año, Carmen Ximénez y Javier Revuelta, como coordinadores del Máster de Metodología de las Ciencias del Comportamiento y de la Salud, deben decidir cuántos Trabajos de Fin de Máster (TFM) de investigación y de prácticas van a ofertar.
Para ello, necesitan conocer la proporción esperada de estudiantes que optarán por un TFM de investigación. Denotemos esta proporción como \(\theta\), un valor desconocido que podemos estimar a partir de los datos del curso anterior.
En esta situación, el método de los momentos es la media aritmética de la variable TFM.
Como \(\hat\theta=0.6\), el 60% de los alumnos escogieron un TFM de investigación el año pasado.
\[ \Pr\left(\text{TFM}_i=1\right)\sim\mbox{Bernoulli}\left(\theta\right) \]
R, la ecuación anterior se representa con dbinom, indicando size=1.\[ \Pr\left(D \mid \theta\right) = \prod^N_{i=1}\mbox{Bernoulli}\left(\text{TFM}_i\mid\theta\right) \]
6.871948e-10 representa la verosimilitud de los datos dado un valor de \(\theta = 0.2\). ¿Qué tal si probamos ahora \(\theta = 0.3\)?3.063652e-08 es mayor que 6.871948e-10, \(\theta=0.3\) explica mejor los datos que \(\theta=0.2\).¿Qué pasa si repetimos este mismo proceso para todos los valores de \(\theta\) entre 0 y 1?
# Valores de theta de 0 a 1
theta <- seq(0, 1, 0.01)
# Vector vacío para guardar la verosimilitud con cada theta
likelihood <- vector(length = length(theta))
# Verosimilitud para cada theta
for(i in 1:length(theta)){
likelihood[i] <- prod(dbinom(TFM, size = 1, prob = theta[i]))
}
# Gráfico de densidad
plot(theta, likelihood)\[ \Pr\left(\theta \mid D\right) \propto \Pr\left(D \mid \theta\right) \times \Pr(\theta) \]
Note
El símbolo \(\propto\) significa “proporcional a”, y no “igual a”. Esto se debe a que la distribución posterior es una distribución de probabilidad, es decir, la suma de sus valores es 1. Sin embargo, esto no se cumple con el producto de la verosimilitud y la distribución a priori.
\[ \Pr\left(\theta \mid D\right) = \frac{\Pr\left(D \mid \theta\right) \times \Pr(\theta)}{\Pr(D)} \]
\[ \Pr\left(\theta \mid D\right) = \frac{\Pr\left(D \mid \theta\right) \times \Pr(\theta)}{\Pr(D)} \]
\[ \Pr\left(\theta \mid D\right) = \frac{\Pr\left(D \mid \theta\right) \times \Pr(\theta)}{\Pr(D)} \]
| Selección de TFM de prácticas e investigación | |||
|---|---|---|---|
| Curso | Investigación | Prácticas | \(\hat{\theta}\) |
| 2018/2019 | 10 | 14 | 0.417 |
| 2019/2020 | 5 | 15 | 0.250 |
| 2020/2021 | 17 | 13 | 0.567 |
| 2021/2022 | 13 | 15 | 0.464 |
| 2022/2023 | 8 | 18 | 0.308 |
\[ \theta \sim \mathcal{N}\left(\mu_\theta = 0.4,\, \sigma_\theta = 0.12\right) \]
R muestra cómo surge toda la magia que crea la famosa distribución posterior.\[ \Pr\left(\theta \mid D\right) = \frac{\Pr\left(D \mid \theta\right) \times \Pr(\theta)}{\Pr(D)}=\frac{\text{Verosimilitud }\times\text{ Prior}}{\text{Verosimilitud marginal}} \]
# Distribución a priori de theta
prior <- dnorm(theta, mean = 0.4, sd = 0.12)
# Distribución posterior (proporcional)
posterior_prop <- likelihood * prior
# Verosimilitud marginal
marginal <- sum(likelihood * prior)
# Distribución posterior escalada a [0,1]
posterior <- posterior_prop/marginal
plot(theta, posterior) # Gráfico de densidad posterior\[ \Pr\left(\theta \mid D\right) = \frac{\Pr\left(D \mid \theta\right) \times \color{darkorange}{\Pr(\theta)}}{\Pr(D)}=\frac{\text{Verosimilitud }\times\color{darkorange}{\text{ Prior}}}{\text{Verosimilitud marginal}} \]
# Distribución a priori de theta
prior <- dnorm(theta, mean = 0.4, sd = 0.12)
# Distribución posterior (proporcional)
posterior_prop <- likelihood * prior
# Verosimilitud marginal
marginal <- sum(likelihood * prior)
# Distribución posterior escalada a [0,1]
posterior <- posterior_prop/marginal
plot(theta, posterior) # Gráfico de densidad posterior\[ \Pr\left(\theta \mid D\right) = \frac{\color{darkorange}{\Pr\left(D \mid \theta\right) \times \Pr(\theta)}}{\Pr(D)}=\frac{\color{darkorange}{\text{Verosimilitud }\times\text{ Prior}}}{\text{Verosimilitud marginal}} \]
# Distribución a priori de theta
prior <- dnorm(theta, mean = 0.4, sd = 0.12)
# Distribución posterior (proporcional)
posterior_prop <- likelihood * prior
# Verosimilitud marginal
marginal <- sum(likelihood * prior)
# Distribución posterior escalada a [0,1]
posterior <- posterior_prop/marginal
plot(theta, posterior) # Gráfico de densidad posterior\[ \Pr\left(\theta \mid D\right) = \frac{\Pr\left(D \mid \theta\right) \times \Pr(\theta)}{\color{royalblue}{\Pr(D)}}=\frac{\text{Verosimilitud }\times\text{ Prior}}{\color{royalblue}{\text{Verosimilitud marginal}}} \]
# Distribución a priori de theta
prior <- dnorm(theta, mean = 0.4, sd = 0.12)
# Distribución posterior (proporcional)
posterior_prop <- likelihood * prior
# Verosimilitud marginal
marginal <- sum(likelihood * prior)
# Distribución posterior escalada a [0,1]
posterior <- posterior_prop/marginal
plot(theta, posterior) # Gráfico de densidad posterior\[ \color{red}{\Pr\left(\theta \mid D\right)} = \frac{\Pr\left(D \mid \theta\right) \times \Pr(\theta)}{\Pr(D)}=\frac{\text{Verosimilitud }\times\text{ Prior}}{\text{Verosimilitud marginal}} \]
# Distribución a priori de theta
prior <- dnorm(theta, mean = 0.4, sd = 0.12)
# Distribución posterior (proporcional)
posterior_prop <- likelihood * prior
# Verosimilitud marginal
marginal <- sum(likelihood * prior)
# Distribución posterior escalada a [0,1]
posterior <- posterior_prop/marginal
plot(theta, posterior) # Gráfico de densidad posterior\[ \theta \sim \mathcal{N}\left(\mu_\theta = \color{royalblue}{0.5},\, \sigma_\theta = \color{royalblue}{0.1}\right) \]
\[ \theta \sim \mathcal{N}\left(\mu_\theta = \color{royalblue}{0.5},\, \sigma_\theta = \color{royalblue}{0.1}\right) \]
\[ \theta \sim \mathcal{N}\left(\mu_\theta = \color{royalblue}{0.5},\, \sigma_\theta = \color{royalblue}{1}\right) \]
\[ \theta \sim \mathcal{N}\left(\mu_\theta = \color{royalblue}{0.5},\, \sigma_\theta = \color{royalblue}{5}\right) \]
\[ \theta \sim \mathcal{N}\left(\mu_\theta = \color{royalblue}{0.5},\, \sigma_\theta = \color{royalblue}{5}\right) \]
\[ \theta \sim \mathcal{N}\left(\mu_\theta = \color{royalblue}{0.5},\, \sigma_\theta = \color{royalblue}{5}\right) \]
\[ \theta\sim\mbox{Uniform}\left(a = 0,\, b = 1\right) \]
[1] 1 1 1 1 1 1 1 1 1 1 1
\[ \Pr\left(\theta \mid D\right) \propto \Pr\left(D \mid \theta\right) \times \color{red}{\Pr(\theta)} \]
[1] 1 1 1 1 1 1 1 1 1 1 1
\[ \Pr\left(\theta \mid D\right) \propto \Pr\left(D \mid \theta\right) \times \color{red}{\boldsymbol{1}}\\ \]
[1] 1 1 1 1 1 1 1 1 1 1 1
\[ \Pr\left(\theta \mid D\right) \color{red}= \Pr\left(D \mid \theta\right) \]
Universo frecuentista, multiverso bayesiano
La estadística frecuentista sólo es un caso particular dentro de la estadística bayesiana.
\[ \color{seagreen}{Y_i} = \color{#0197FD}{\beta_0 + \beta_1\cdot X_i} + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\color{#7E57C2}{\mathcal{N}}\left(\color{#9a2515}\mu =0,\, \color{#9a2515}{\sigma} = \color{#0197FD}{\sigma_\varepsilon}\right) \]
\[ \color{seagreen}{Y_i} \sim \color{#7E57C2}{\mathcal{N}}\left(\color{#9a2515}\mu = \color{#0197FD}{\beta_0 + \beta_1\cdot X_i},\,\color{#9a2515}{\sigma} = \color{#0197FD}{\sigma_\varepsilon} \right) \]
\[ Y_i = \beta_0 + \beta_1\cdot X_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \sigma_\varepsilon\right) \]
\[ Y_i = \nu + \lambda \cdot \eta_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta\right) \]
\[ Y_i = \nu + \lambda \cdot \eta_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta\right) \]
\[ \begin{aligned} & Y_{i1} = \nu_1 + \lambda_1 \cdot \eta_i + \varepsilon_{i1},\quad\text{donde}\quad\varepsilon_{i1}\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta_1\right) \end{aligned} \]
\[ Y_i = \nu + \lambda \cdot \eta_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta\right) \]
\[ \begin{aligned} & Y_{i1} = \nu_1 + \lambda_1 \cdot \eta_i + \varepsilon_{i1},\quad\text{donde}\quad\varepsilon_{i1}\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta_1\right)\\ & Y_{i2} = \nu_2 + \lambda_2 \cdot \eta_i + \varepsilon_{i2},\quad\text{donde}\quad\varepsilon_{i2}\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta_2\right) \end{aligned} \]
\[ Y_i = \nu + \lambda \cdot \eta_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta\right) \]
\[ \begin{aligned} & Y_{i1} = \nu_1 + \lambda_1 \cdot \eta_i + \varepsilon_{i1},\quad\text{donde}\quad\varepsilon_{i1}\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta_1\right)\\ & Y_{i2} = \nu_2 + \lambda_2 \cdot \eta_i + \varepsilon_{i2},\quad\text{donde}\quad\varepsilon_{i2}\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta_2\right)\\ & Y_{i3} = \nu_3 + \lambda_3 \cdot \eta_i + \varepsilon_{i3},\quad\text{donde}\quad\varepsilon_{i3}\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta_3\right)\\ \end{aligned} \]
\[ Y_i = \nu + \lambda \cdot \eta_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta\right) \]
\[ \begin{aligned} & Y_{i1} = \nu_1 + \lambda_1 \cdot \color{royalblue}{\eta_i} + \varepsilon_{i1},\quad\text{donde}\quad\varepsilon_{i1}\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta_1\right)\\ & Y_{i2} = \nu_2 + \lambda_2 \cdot \color{royalblue}{\eta_i} + \varepsilon_{i2},\quad\text{donde}\quad\varepsilon_{i2}\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta_2\right)\\ & Y_{i3} = \nu_3 + \lambda_3 \cdot \color{royalblue}{\eta_i} + \varepsilon_{i3},\quad\text{donde}\quad\varepsilon_{i3}\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta_3\right)\\ \end{aligned} \]
\[ Y_i = \nu + \lambda \cdot \eta_i + \varepsilon_i,\quad\text{donde}\quad\varepsilon_i\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta\right) \]
\[ Y_{i\color{red}j} = \nu_\color{red}j + \lambda_\color{red}j \cdot \eta_i + \varepsilon_{i\color{red}j},\quad\text{donde}\quad\varepsilon_{i\color{red}j}\sim\mathcal{N}\left(\mu =0,\, \sigma = \theta_\color{red}j\right) \]
\[ \begin{aligned} & Y_{ij} = \nu_j + \color{royalblue}{\lambda_{j1} \cdot \eta_{i1}} + \varepsilon_{ij}\\ \end{aligned} \]
\[ \begin{aligned} & Y_{ij} = \nu_j + \color{royalblue}{\lambda_{j1} \cdot \eta_{i1}} + \varepsilon_{ij}\\ & Y_{ij} = \nu_j + \color{royalblue}{\lambda_{j1} \cdot \color{royalblue}\eta_{i1}} + \color{firebrick}{\lambda_{j2} \cdot \eta_{i2}} + \varepsilon_{ij}\\ \end{aligned} \]
\[ \begin{aligned} & Y_{ij} = \nu_j + \color{royalblue}{\lambda_{j1} \cdot \eta_{i1}} + \varepsilon_{ij}\\ & Y_{ij} = \nu_j + \color{royalblue}{\lambda_{j1} \cdot \eta_{i1}} + \color{firebrick}{\lambda_{j2} \cdot \eta_{i2}} + \varepsilon_{ij}\\ & Y_{ij} = \nu_j + \color{royalblue}{\lambda_{j1} \cdot \eta_{i1}} + \color{firebrick}{\lambda_{j2} \cdot \eta_{i2}} + \color{seagreen}{\lambda_{j3} \cdot \eta_{i3}} + \varepsilon_{ij}\\ \end{aligned} \]
\[ \begin{aligned} & Y_{ij} = \nu_j + \color{royalblue}{\lambda_{j1} \cdot \eta_{i1}} + \varepsilon_{ij}\\ & Y_{ij} = \nu_j + \color{royalblue}{\lambda_{j1} \cdot \eta_{i1}} + \color{firebrick}{\lambda_{j2} \cdot \eta_{i2}} + \varepsilon_{ij}\\ & Y_{ij} = \nu_j + \color{royalblue}{\lambda_{j1} \cdot \eta_{i1}} + \color{firebrick}{\lambda_{j2} \cdot \eta_{i2}} + \color{seagreen}{\lambda_{j3} \cdot \eta_{i3}} + \varepsilon_{ij}\\ & Y_{ij} = \nu_j + \color{royalblue}{\lambda_{j1} \cdot \eta_{i1}} + \color{firebrick}{\lambda_{j2} \cdot \eta_{i2}} + \color{seagreen}{\lambda_{j3} \cdot \eta_{i3}} + \dots + \varepsilon_{ij}\\ \end{aligned} \]
\[ Y_{ij} = \nu_j + \sum^ \color{red}M_{ \color{red}m=1} \lambda_{j \color{red}m} \cdot \eta_{i \color{red}m} + \varepsilon_{ij} \]
\[ \underbrace{ \begin{bmatrix} Y_{i1}\\ Y_{i2}\\ Y_{i3}\\ \end{bmatrix}}_{\textstyle{\mathbf{y}_i}} = \underbrace{ \begin{bmatrix} \nu_{1}\\ \nu_{2}\\ \nu_{3}\\ \end{bmatrix}}_{\textstyle{\boldsymbol{\nu}}} + \underbrace{ \begin{bmatrix} \lambda_{11} & \lambda_{12}\\ \lambda_{21} & \lambda_{22}\\ \lambda_{31} & \lambda_{32}\\ \end{bmatrix}}_{\textstyle{\mathbf{\Lambda}}} \times \underbrace{ \begin{bmatrix} \eta_{i1}\\ \eta_{i2}\\ \end{bmatrix}}_{\textstyle{\boldsymbol{\eta}_i}} + \underbrace{ \begin{bmatrix} \varepsilon_{1}\\ \varepsilon_{2}\\ \varepsilon_{3}\\ \end{bmatrix}}_{\textstyle{\boldsymbol{\varepsilon}_i}} \]
\[ \color{seagreen}{Y_{ij}} \sim \color{#7E57C2}{\mathcal{N}}\left(\color{#9a2515}\mu = \color{#0197FD}{\nu_j + \sum^M_{m=1}\lambda_{jm} \cdot \eta_{im}},\,\color{#9a2515}\sigma = \color{#0197FD}{\theta}\right) \]
\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]
\[ Y_{ij} \sim \mathcal{N}\left(\mu = \nu_j + \sum^M_{m=1}\lambda_{jm} \cdot \color{firebrick}{\eta_{im}},\,\sigma = \theta\right) \]
\[ \mathbf{y}_i \sim \mathcal{MVN}\left(\boldsymbol\mu = \boldsymbol\nu + \mathbf{\Lambda} \cdot \color{firebrick}{\boldsymbol{\alpha}},\,\mathbf{\Sigma} = \mathbf{\Lambda\cdot\color{firebrick}{\mathbf\Psi}\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}\right) \]
RModelo sin estructura de medias
Modelo con estructura de medias
R\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]
\[ \boldsymbol{\nu} = \left[-1,\,0,\,1,\,-1,\,0,\,1\right] \]
R\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]
\[ \mathbf{\Lambda}^\prime=\begin{bmatrix} 1.40 & 1.75 & 2.10 & 0 & 0 & 0\\ 0 & 0 & 0 & 1.40 & 1.75 & 2.10 \\ \end{bmatrix} \]
R\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]
\[ \boldsymbol{\alpha} = \left[0,\,0\right] \]
R\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]
\[ \mathbf{\Psi} = \begin{bmatrix} 1 & 0.5\\ 0.5 & 1 \end{bmatrix} \]
R\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]
\[ \mathbf{\Theta}=\begin{bmatrix} \theta^2_{11} & 0 & 0 & 0 & 0 & 0\\ 0 & \theta^2_{22} & 0 & 0 & 0 & 0\\ 0 & 0 & \theta^2_{33} & 0 & 0 & 0\\ 0 & 0 & 0 & \theta^2_{44} & 0 & 0\\ 0 & 0 & 0 & 0 & \theta^2_{55} & 0\\ 0 & 0 & 0 & 0 & 0 & \theta^2_{66}\\ \end{bmatrix} \]
R\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]
\[ \mathbf{\Theta}=\begin{bmatrix} 2.04 & 0 & 0 & 0 & 0 & 0\\ 0 & 3.1875 & 0 & 0 & 0 & 0\\ 0 & 0 & 4.59 & 0 & 0 & 0\\ 0 & 0 & 0 & 2.04 & 0 & 0\\ 0 & 0 & 0 & 0 & 3.1875 & 0\\ 0 & 0 & 0 & 0 & 0 & 4.59\\ \end{bmatrix} \]
R\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]
# Matriz de covarianzas únicas poblacional
(Theta <- diag(c(2.04, 3.1875, 4.59, 2.04, 3.1875, 4.59), 6, 6)) [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 2.04 0.0000 0.00 0.00 0.0000 0.00
[2,] 0.00 3.1875 0.00 0.00 0.0000 0.00
[3,] 0.00 0.0000 4.59 0.00 0.0000 0.00
[4,] 0.00 0.0000 0.00 2.04 0.0000 0.00
[5,] 0.00 0.0000 0.00 0.00 3.1875 0.00
[6,] 0.00 0.0000 0.00 0.00 0.0000 4.59
R\[ \color{seagreen}{\mathbf{y}_i} \sim \color{#7E57C2}{\mathcal{MVN}}\left(\color{#9a2515}{\boldsymbol\mu} = \color{#0197FD}{\boldsymbol\nu + \mathbf{\Lambda} \cdot \boldsymbol{\alpha}},\,\color{#9a2515}{\mathbf{\Sigma}} = \color{#0197FD}{\mathbf{\Lambda\cdot\mathbf\Psi\cdot\mathbf{\Lambda}^\prime + \mathbf{\Theta}}}\right) \]
{blavaan}{blavaan} es idéntica a la de {lavaan}.{blavaan} para los parámetros del modelo?| Parámetro | priors | Parámetro | priors |
|---|---|---|---|
| \(\nu\) | normal(0,32) | \(\theta\) | gamma(1,.5)[sd] |
| \(\alpha\) | normal(0,10) | \(\psi\) | gamma(1,.5)[sd] |
| \(\lambda\) | normal(0,10) | \(\rho\) | beta(1,1) |
{blavaan} tiene un atajo para realizar prior predictive checks.{bayesplot}{blavaan} no se puede cambiar la distribución a priori, sólo sus parámetros.# Nueva distribución a priori
(new_priors <- dpriors(nu = "normal(0,15)",
lambda = "normal(0,5)",
rho = "beta(2,2)",
psi = "gamma(1,.2)"))
# Muestreamos de la nueva distribución a priori
new_prior_pred <- bcfa(model.2f, data = Y, std.lv = TRUE,
meanstructure = TRUE, test = "none",
sample = 2e4, prisamp = TRUE,
dp = new_priors, # Added new priors
bcontrol = list(cores = 3)) Estimaremos los modelos de uno y dos factores con las nuevas priors.
Utilizaremos 3 cadenas de Markov con 1500 iteraciones, descartando las primeras 500 iteraciones.
bcontrol permite paralelizar las cadenas de Markov.
# Modelo de un factor
bcfa.1f.fit <- bcfa(model.1f, data = Y, burnin = 500,
sample = 1000, meanstructure = TRUE,
std.lv = TRUE, bcontrol = list(cores = 3))
# Modelo de dos factores
bcfa.2f.fit <- bcfa(model.2f, data = Y, burnin = 500,
sample = 1000, meanstructure = TRUE,
std.lv = TRUE, bcontrol = list(cores = 3))Índices de ajuste incremental en {blavaan}
Los índices de ajuste incremental bayesianos, BCFI y BTLI, comparan el ajuste del modelo propuesto con el de un modelo nulo. Sin embargo, {blavaan} no estima ese modelo nulo por defecto, como sí hace {lavaan}. Por ello, debemos especificar manualmente el modelo nulo y estimarlo.
EAP Median MAP SD lower upper
BRMSEA 0.313 0.313 0.312 0.007 0.301 0.327
BGammaHat 0.873 0.874 0.874 0.005 0.864 0.882
adjBGammaHat 0.228 0.230 0.235 0.030 0.169 0.282
BMc 0.804 0.805 0.806 0.008 0.789 0.818
BCFI 0.714 0.715 0.716 0.013 0.688 0.738
BTLI 0.052 0.056 0.059 0.044 -0.035 0.131
BNFI 0.714 0.715 0.716 0.013 0.688 0.737
EAP Median MAP SD lower upper
BRMSEA 0.046 0.050 0.057 0.026 0.000 0.085
BGammaHat 0.993 0.993 1.000 0.006 0.981 1.000
adjBGammaHat 0.975 0.978 0.998 0.021 0.936 1.000
BMc 0.989 0.990 0.999 0.009 0.972 1.000
BCFI 0.985 0.987 0.999 0.012 0.962 1.000
BTLI 0.974 0.976 0.983 0.024 0.924 1.017
BNFI 0.969 0.970 0.974 0.013 0.943 0.992
{blavaan} incluye dos métodos para comparar modelos:{loo}.{blavaan} no recomienda su uso.
Un factor
|
Dos factores
|
|||
|---|---|---|---|---|
| Estimación | Error típico | Estimación | Error típico | |
| elpd_loo | -4048.296 | 34.544 | -3984.390 | 32.102 |
| p_loo | 22.568 | 1.663 | 18.993 | 1.202 |
| looic | 8096.592 | 69.087 | 7968.780 | 64.204 |
Model2 es el modelo que hemos indicado previamente en object 2
Latent Variables:
Estimate Std.lv Std.all
F1 =~
V1 1.452 1.452 0.711
V2 2.091 2.091 0.759
V3 2.031 2.031 0.673
F2 =~
V4 1.566 1.566 0.722
V5 1.671 1.671 0.708
V6 1.905 1.905 0.667
Covariances:
Estimate Std.lv Std.all
F1 ~~
F2 0.469 0.469 0.469
Intercepts:
Estimate Std.lv Std.all
.V1 -1.086 -1.086 -0.532
.V2 -0.013 -0.013 -0.005
.V3 1.029 1.029 0.341
.V4 -0.918 -0.918 -0.423
.V5 0.041 0.041 0.017
.V6 0.986 0.986 0.345
F1 0.000 0.000 0.000
F2 0.000 0.000 0.000
Variances:
Estimate Std.lv Std.all
.V1 2.060 2.060 0.494
.V2 3.218 3.218 0.424
.V3 4.974 4.974 0.547
.V4 2.247 2.247 0.478
.V5 2.779 2.779 0.499
.V6 4.538 4.538 0.556
F1 1.000 1.000 1.000
F2 1.000 1.000 1.000
{blavaan} podemos utilizar la función sampleData para generar datos desde el modelo.# Simulamos 100 bases de datos
yrep <- sampleData(bcfa.2f.fit, nrep = 100, simplify = TRUE)
# Guardamos los resultados para cada item
posterior_plots_list <- vector(mode = "list", length = ncol(Y))
for(i in 1:ncol(Y)){
posterior_plots_list[[i]] <- ppc_dens_overlay(
y = Y[,i], yrep = t(sapply(yrep, function(x) x[,i]))
) +
labs(title = paste("Item", i))
}
# Gráfico comparando y e yrep
cowplot::plot_grid(plotlist = posterior_plots_list){blavaan} pueden realizarse con la función ppmcppmc puede utilizarse para calcular el estadístico \(\chi^2\) y SRMR en los datos observados y simulados en cada iteración.plot e hist comparan gráficamente los datos observados y simulados.ppmc con cualquier función que pueda utilizarse para un modelo de {lavaan}.compRelSEM de {SemTools} para obtener la distribución posterior de \(\omega\).# Función para obtener la fiabilidad en cada iteración
bcfa.reliability <- function(fit){ semTools::compRelSEM(fit) }
# Fiabilidad en el modelo unidimensional
ppmc_omega_1f <- ppmc(object = blavaan.1f.fit,
discFUN = bcfa.reliability)
# Fiabilidad en el modelo de dos factores
ppmc_omega_2f <- ppmc(object = blavaan.2f.fit,
discFUN = bcfa.reliability){cmdstanr}R, podemos usar {Rstan} o {cmdstanr}.{cmdstanr}:
{cmdstanr} utiliza la última versión de Stan.{posterior}, {loo}, {bayesplot}…{cmdstanr}.Rstudio.R como una lista.// Espacio para funciones
functions {
}
// Espacio para los datos
data {
}
// Espacio para transformar los datos
transformed data {
}
// Espacio para los parámetros
parameters {
}
// Espacio para los parámetros transformados
transformed parameters {
}
// Espacio para el modelo
model {
}
// Espacio para las cantidades generadas
generated quantities{
}\[ \begin{aligned} &\boldsymbol\mu\sim\mbox{Normal}\left(\mu = 0,\,\, \sigma = 10\right)\\ &\mathbf\Lambda\sim\mbox{Normal}\left(\mu = 0,\,\, \sigma = 5\right)\\ &\boldsymbol\theta\sim\mbox{Gamma}\left(\alpha = 1,\,\, \beta = 0.5\right)\\ &L_\Psi\sim\mbox{LkjCholesky}\left(\eta=1\right)\\ \end{aligned} \]
\[ \mathbf{y}_i\sim\mathcal{MVN}\left(\boldsymbol\mu,\,\,\mathbf{\Sigma}\right) \]
{blavaan} soluciona este problema de signos en las cantidades generadas.{cmdstanr}{cmdstanr} creará un programa ejecutable en vuestro ordenador.R.BCFA es una especie de lista con todas las funciones que pueden utilizarse dentro.{cmdstanr}BCFA$sample() es la función para estimar el modelo bayesiano.{cmdstanr}BCFA$sample() cuenta con muchos argumentos, pero no son similares a los de {blavaan}BCFA_1f <- BCFA$sample(
data = sdata.1f, # Stan data
chains = 4, # Number of chains
parallel_chains = 4, # Number of parallel chains
iter_warmup = 500, # Adaptation iterations
iter_sampling = 1500, # Sampled iterations
refresh = 500, # Progress bar at 500 iterations
init = 0) # All starting values = 0BCFA-1f será una lista que contenga los resultados y funciones para almacenarlos y resumirlos.{cmdstanr}BCFA_1f$summary() genera tablas resumen de los parámetros que indiquemos.| variable | mean | median | sd | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|---|
| Lambda[1,1] | 0.72 | 0.72 | 0.13 | 0.13 | 0.52 | 0.93 | 1 | 8877.06 | 5136.89 |
| Lambda[2,1] | 1.17 | 1.17 | 0.17 | 0.17 | 0.90 | 1.45 | 1 | 6893.36 | 5138.69 |
| Lambda[3,1] | 1.30 | 1.29 | 0.20 | 0.20 | 0.97 | 1.63 | 1 | 6140.56 | 4975.74 |
| Lambda[4,1] | 1.53 | 1.53 | 0.13 | 0.13 | 1.32 | 1.74 | 1 | 6879.88 | 4256.78 |
| Lambda[5,1] | 1.99 | 1.99 | 0.16 | 0.16 | 1.72 | 2.26 | 1 | 7426.89 | 4746.85 |
| Lambda[6,1] | 1.89 | 1.89 | 0.18 | 0.18 | 1.58 | 2.19 | 1 | 7576.66 | 4991.18 |
{cmdstanr}{cmdstanr}BCFA_2f <- BCFA$sample(
data = sdata.2f, # Stan data
chains = 4, # Number of chains
parallel_chains = 4, # Number of parallel chains
iter_warmup = 500, # Adaptation iterations
iter_sampling = 1500, # Sampled iterations
refresh = 500, # Progress bar
seed = 2025, # Reproducible results
init = 0) # All starting values = 0{cmdstanr}| variable | mean | median | sd | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|---|
| Lambda[1,1] | 1.09 | 1.09 | 0.12 | 0.12 | 0.90 | 1.29 | 1 | 6068.92 | 4459.71 |
| Lambda[2,1] | 1.71 | 1.71 | 0.17 | 0.17 | 1.44 | 1.98 | 1 | 4970.63 | 4170.68 |
| Lambda[3,1] | 2.20 | 2.20 | 0.18 | 0.19 | 1.90 | 2.50 | 1 | 4619.45 | 4105.59 |
| Lambda[4,2] | 1.57 | 1.57 | 0.12 | 0.12 | 1.37 | 1.78 | 1 | 5773.14 | 4878.45 |
| Lambda[5,2] | 2.05 | 2.05 | 0.16 | 0.16 | 1.78 | 2.31 | 1 | 5678.16 | 4705.00 |
| Lambda[6,2] | 2.02 | 2.01 | 0.18 | 0.18 | 1.73 | 2.31 | 1 | 5434.57 | 4405.27 |
$draws(){loo} elpd_diff se_diff
model2 0.0 0.0
model1 -55.2 12.0
Ojalá esto fuese más sencillo
La comparación de modelos en estadística bayesiana ha sido (y sigue siendo) un tema de fervoroso debate entre dos disciplinas: los psicólogos matemáticos y los estadísticos bayesianos.
\[ \Pr\left(\theta \mid D\right) = \frac{\Pr\left(D \mid \theta\right) \times \Pr(\theta)}{\color{royalblue}{\Pr(D)}} \]
\[ \mbox{BF}_{12}=\frac{\color{royalblue}{\Pr(D\mid\mathcal{M}_1)}}{\color{royalblue}{\Pr(D\mid\mathcal{M}_2)}} \]
\[ \mbox{BF}_{12}=\frac{\Pr(D\mid\mathcal{M}_1)}{\Pr(D\mid\mathcal{M}_2)} \]
for (i in 1:n) {
# Ajuste sin la i-ésima observación
fit_loo <- lm(y[-i] ~ x[-i])
# Predicción para la observación i
y_pred_i <- predict(fit_loo, newdata = data.frame(x = x[i]))
# Error de estimación sin la i-ésima obs
sigma_i <- summary(fit_loo)$sigma
# Cálculo del log-densidad predictiva de y[i]
log_pdens[i] <- dnorm(y[i], mean = y_pred_i, sd = sigma_i,
log = TRUE)
}
# ELPD: suma de las log-densidades
ELPD_loo <- sum(log_pdens)ELPD son las siglas de Expected log pointwise predictive density, y su valor refleja la calidad predictiva del modelo.
La comparación entre ELPDs permite saber qué modelo predice mejor los datos fuera de la muestra.
\[ \Delta\mbox{ELPD}_{\mathcal{M}_1,\,\mathcal{M}_2}=\mbox{ELPD}_{\mathcal{M}_1} - \mbox{ELPD}_{\mathcal{M}_2} \]
{loo} estima \(\Delta\mbox{ELPD}_{\mathcal{M}_1,\,\mathcal{M}_2}\) y su error típico.\[ \mbox{PsBF}_{12}=\exp\left(\Delta\mbox{ELPD}_{\mathcal{M}_1,\,\mathcal{M}_2}\right) \]
{blavaan} incluye una estimación de la verosimilitud marginal y del factor de bayes basada en la aproximación de Laplace-Metrópolis.Los datos observados son \(1.16\times 10^{27}\) veces más probables con el modelo de dos factores que con el modelo de un factor.
Finalmente, el Pseudo-Factor de Bayes se obtiene exponenciando \(\Delta\mbox{ELPD}\).